#include <iostream>
#include <filesystem>

#include "gdal/ogr_geometry.h"
#include "gdal/gdal_priv.h"
#include "gdal/gdal_alg.h"

#include"Common.cuh"
#include "../CuLibDGGSKernel/ICuDGGSKernel.cuh"
#include "LibGridTIFF/GridTIFF.hpp"

using namespace NNU::DGGS;
using namespace NNU::DGGS::Base;
using namespace NNU::DGGS::Kernel;
namespace fs = std::filesystem;

#define RasterIndex std::string

/// GPU相关参数
int pGridDimX=0;
int pBlockDimX=0;

/// 格网相关参数
DGGSGridType pGridType=DGGSGridType::Unknown;
DGGSElementType pElementType=DGGSElementType::Unknown;

/// CUDGGS的内核
ICuDGGSKernel* pCuKernel= nullptr;

///raster映射多个GridTIFF,用于Raster2Grid
std::multimap<RasterIndex,DGGSGridTileCode>*raster2TileMap= nullptr;

///合并相关参数
bool minMaxFlag=false;
double minValue=0.0;
double maxValue=0.0;
double mergeBound[4]={-90,90,180,-180};
int mergeBandCount=0;
double mergeDiffer=0.0;
GDALDataType mergeDataType=GDT_Unknown;

/*********************************************读写GeoTIFF****************************************************/
#pragma region  GeoTIFF读写声明
/// <summary>
///写入GeoTIFF文件
/// </summary>
/// <param name="fileName">文件名</param>
/// <param name="geoTransform">地理变换</param>
/// <param name="xSize">x轴大小</param>
/// <param name="ySize">y轴大小</param>
/// <param name="dataType">数据类型</param>
/// <param name="pDataList">数据列表</param>
/// <param name="noDataValue">无值</param>
/// <param name="bandCount">波段数</param>
/// <returns></returns>
template<class T>
bool saveRaster(const char* fileName, double* geoTransform, int xSize, int ySize,GDALDataType dataType, T** pDataList,double noDataValue,int bandCount);

/// <summary>
/// 裁剪跨180的菱形
/// </summary>
/// <param name="fileName">文件名</param>
/// <returns></returns>
bool clip(const char *fileName);

/// <summary>
/// 读取GeoTIFF文件
/// </summary>
/// <param name="fileName">文件名</param>
/// <param name="geoTransform">地理变换</param>
/// <param name="xSize">x</param>
/// <param name="ySize">y</param>
/// <param name="dataType">数据类型</param>
/// <param name="pDataList">数据列表</param>
/// <param name="noDataValue">无值</param>
/// <param name="bandCount">波段数</param>
/// <returns></returns>
bool loadRaster(const char* fileName, double* &geoTransform, int &xSize, int &ySize,GDALDataType &dataType, void** &pDataList,double &noDataValue,int &bandCount);

#pragma endregion
/*********************************************Raster2Grid***************************************************/
#pragma region  Raster2Grid声明

/// <summary>
/// 创建Map
/// </summary>
/// <param name="srcFileNames">源文件名</param>
/// <param name="fileCount">文件数量</param>
/// <param name="level">等级</param>
/// <param name="originLevel">原点等级</param>
/// <returns></returns>
DGGSDataType createMap(std::string* srcFileNames,int fileCount,int &level,int originLevel);

/// <summary>
/// 获取点对应的栅格值
/// </summary>
/// <pram name="coord">经纬度</param>
/// <param name="trans">地理变换参数</param>
/// <param name="xSize">栅格宽度</param>
/// <param name="ySize">栅格高度</param>
/// <param name="srcValueList">栅格数据</param>
/// <returns>栅格值</returns>
template<class T>
__device__ T getValueFromCoord(GeoCoord coord, const double* trans, int xSize, int ySize, T* srcValueList,double noValue);

/// <summary>
/// 从栅格中获取GridTIFF的属性值
/// </summary>
/// <param name="srcValueList">栅格数据</param>
/// <param name="trans">地理变换参数</param>
/// <param name="xSize">栅格宽度</param>
/// <param name="ySize">栅格高度</param>
/// <param name="tagValueList">GridTIFF的属性值</param>
/// <param name="tagTrans">GridTIFF的变换参数</param>
/// <param name="gridSize">GridTIFF的栅格尺寸</param>
/// <param name="gridType">GridTIFF的栅格类型</param>
/// <returns></returns>
template<class T>
__global__ void setGridTIFFValuesFromRaster(T* srcValueList, double* trans, int xSize, int ySize, T* tagValueList, GridTransform tagTrans, int gridSize, ICuDGGSKernel* pKernel, double noValue);

/// <summary>
///单张Raster文件为GridTIFF赋值
/// <summary>
/// <param name="srcFileName">源文件名</param>
/// <param name="originLevel">原点等级</param>
/// <param name="tagFilePath">目标文件名</param>
/// <returns></returns>
template<class T>
void singleRaster2Grid(const std::string& srcFileName, const char *exFormat,int originLevel,const std::string& tagFilePath);

/// <summary>
/// 多张Raster文件为GridTIFF赋值
/// </summary>
/// <param name="srcFilePath">源文件路径</param>
/// <param name="tagFilePath">目标文件路径</param>
/// <param name="exFormat">目标文件格式</param>
/// <param name="gridDimX">栅格尺寸</param>
/// <param name="blockDimX">块尺寸</param>
/// <param name="gridType">栅格类型</param>
/// <param name="elementType">元素类型</param>
/// <param name="level">等级</param>
/// <param name="originLevel">原点等级</param>
/// <returns></returns>
void raster2Grid(const char* srcFilePath,  const char* tagFilePath, const char*exFormat,int gridDimX, int blockDimX, DGGSGridType gridType, DGGSElementType elementType, int level, int originLevel);

#pragma endregion
/*********************************************tile相关操作***************************************************/
#pragma region tile与bound相关操作声明

/// <summary>
/// 根据经纬度包围盒获取相关GridTIFF
/// </summary>
/// <param name="bound">边界</param>
/// <param name="differ">差值</param>
/// <param name="xSize">栅格尺寸</param>
/// <param name="ySize">栅格尺寸</param>
/// <param name="xStep">x轴步长</param>
/// <param name="yStep">y轴步长</param>
/// <param name="level">等级</param>
/// <param name="originLevel">原点等级</param>
/// <param name="tiffCodes">瓦片编码</param>
/// <param name="pKernel">CuDGGSKernel</param>
/// <returns></returns>
__global__ void getGridTileCodeFromGridCode(const double* bound, double differ, int xSize, int ySize, int xStep, int yStep, int level, int originLevel, DGGSGridTileCode* tiffCodes, ICuDGGSKernel *pKernel);


/// <summary>
/// 根据经纬度包围盒获取相关GridTIFF
/// </summary>
/// <param name="bound">边界</param>
/// <param name="level">等级</param>
/// <param name="originLevel">原点等级</param>
/// <param name="tileCodes">瓦片编码</param>
/// <param name="tiffCount">瓦片数量</param>
/// <returns></returns>
bool getGridTileFromBound(const double *bound, int level, int originLevel, DGGSGridTileCode *&tileCodes, int &tiffCount);

/// <summary>
/// 要素编码转经纬度
/// </summary>
/// <param name="codes">编码</param>
/// <param name="codeCount">编码数量</param>
/// <param name="coords">经纬度</param>
/// <param name="pKernel">ISEA4HImpl</param>
/// <returns></returns>
__global__ void code2Coord(DGGSElementCode * codes, int codeCount, GeoCoord* coords, ICuDGGSKernel*pKernel);

/// <summary>
/// 根据GridTIFF编码获取空间边界
/// </summary>
/// <param name="gridTiffCode">GridTIFF编码</param>
/// <param name="originLevel">原点等级</param>
/// <returns></returns>
DGGridDiamond getSpaceBoundFromGridTileCode(DGGSGridTileCode gridTiffCode, int originLevel);

#pragma endregion
/*********************************************Grid2Raster***************************************************/
#pragma region Grid2Raster声明

/// <summary>
/// 根据GridTIFF设置Raster值
/// </summary>
/// <param name="srcValueList">源数据</param>
/// <param name="srcTrans">源数据变换</param>
/// <param name="tileSize">块尺寸</param>
/// <param name="tagValueList">目标数据</param>
/// <param name="trans">变换</param>
/// <param name="xSize">栅格尺寸</param>
/// <param name="ySize">栅格尺寸</param>
/// <param name="noValue">无效值</param>
/// <param name="pKernel">ISEA4HImpl</param>
/// <returns></returns>
template<class T>
__global__ void setRasterValuesFromGrid(T* srcValueList, GridTransform srcTrans, int tileSize, T* tagValueList, double* trans, int xSize, int ySize, double noValue,double differTol, ICuDGGSKernel* pKernel);

/// <summary>
/// 单张Raster转化
/// </summary>
/// <param name="srcFilePath">源文件路径</param>
/// <param name="gridTiffName">GridTIFF名称</param>
/// <param name="exFormat">输出格式</param>
/// <param name="tagFilePath">目标文件路径</param>
template<class T>
void singleGridTranslate(const char* srcFilePath,const char* gridTiffName,const char*exFormat, const char* tagFilePath);

/// <summary>
/// 合并Raster值
/// </summary>
/// <param name="srcValueList">源数据</param>
/// <param name="srcXSize">源数据X尺寸</param>
/// <param name="srcYSize">源数据Y尺寸</param>
/// <param name="tagValueList">目标数据</param>
/// <param name="noValue">无效值</param>
template<class T>
__global__ void mergeRasterValues(T* srcValueList, int srcXSize,int srcYSize,T*tagValueList,double noValue);

/// <summary>
/// 合并Raster
/// </summary>
/// <param name="srcFilePath">源文件路径</param>
/// <param name="tagFileName">目标文件名</param>
/// <param name="exFormat">输出格式</param>
template<class T>
void mergeRaster(const char* srcFilePath,const char* tagFileName,const char*exFormat);

/// <summary>
/// 单组GridTIFF转化
/// </summary>
/// <param name="srcFilePath">源文件路径</param>
/// <param name="tagFilePath">目标文件路径</param>
/// <param name="exFormat">输出格式</param>
void singleGrid2Raster(const char* srcFilePath, const char* tagFilePath,const char*exFormat);

/// <summary>
/// 多组GridTIFF转化
/// </summary>
/// <param name="srcFilePath">源文件路径</param>
/// <param name="tagFilePath">目标文件路径</param>
/// <param name="exFormat">输出格式</param>
void multGrid2Raster(const char* srcFilePath, const char* tagFilePath,const char*exFormat);

#pragma endregion
/*********************************************读写GeoTIFF***************************************************/
#pragma region GeoTIFF读写实现

template<class T>
bool saveRaster(const char* fileName, double* geoTransform, int xSize, int ySize,GDALDataType dataType, T** pDataList,double noDataValue,int bandCount)
{
    GDALAllRegister();
    const char* pszFormat = "GTiff";
    GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (poDriver == nullptr)
    {
        return false;
    }
    char** papszOptions = nullptr;
    GDALDataset* poDstDS = poDriver->Create(fileName, xSize, ySize, bandCount, dataType, papszOptions);
    poDstDS->SetGeoTransform(geoTransform);
    // 设置投影为WGS84
    OGRSpatialReference srs;
    srs.SetWellKnownGeogCS("WGS84");
    char* wkt = nullptr;
    srs.exportToWkt(&wkt);
    poDstDS->SetProjection(wkt);
    CPLFree(wkt);

    for (int i = 0; i < bandCount; i++)
    {
        poDstDS->GetRasterBand(i+1)->SetNoDataValue(noDataValue);
        if(pDataList[i]!= nullptr)
            poDstDS->GetRasterBand(i+1)->RasterIO(GF_Write, 0, 0, xSize, ySize, pDataList[i], xSize, ySize, dataType, 0, 0);
    }
    GDALClose((GDALDatasetH)poDstDS);

    return true;
}

bool clip(const char *fileName)
{
    // 输入和输出文件路径
    std::string inputFile = fileName;

    //后缀
    std::string exterm;
    size_t pos = inputFile.find_last_of('.');
    if (pos != std::string::npos) {
        exterm= inputFile.substr(pos + 1);
    }

    std::string outputFile1 = inputFile.substr(0, inputFile.find_last_of('.')) + "_1."+exterm;
    std::string outputFile2 = inputFile.substr(0, inputFile.find_last_of('.')) + "_2."+exterm;

    // 初始化GDAL
    GDALAllRegister();

    // 打开输入TIFF文件
    auto *poDataset = (GDALDataset *) GDALOpen(inputFile.c_str(), GA_ReadOnly);
    if (poDataset == nullptr) {
        std::cerr << "Error: Unable to open input file " << inputFile << std::endl;
        return false;
    }

    //获取参数
    double noDataValue = poDataset->GetRasterBand(1)->GetNoDataValue();
    int bandCount= poDataset->GetRasterCount();
    GDALDataType dataType = poDataset->GetRasterBand(1)->GetRasterDataType();
    double adfGeoTransform[6];
    poDataset->GetGeoTransform(adfGeoTransform);
    double minX = adfGeoTransform[0];
    double pixelWidth = adfGeoTransform[1];
    double maxX = minX + (poDataset->GetRasterXSize() * pixelWidth);
    int nXSize = poDataset->GetRasterXSize();
    int nYSize = poDataset->GetRasterYSize();


    double leftMinLon = -180;
    double rightMinLon = 180;
    // 分割的范围

    double rightBound=-180+std::ceil(360.0/pixelWidth)*pixelWidth;;
    int left1 = 0,right1 = std::round((rightBound - minX) / pixelWidth);
    leftMinLon=adfGeoTransform[0];
    rightMinLon=-180.0;

    if(minX<-180)
    {
        double temp =(-180 - minX) / pixelWidth;

        right1= std::round((-180 - minX) / pixelWidth);
        leftMinLon=360.0+adfGeoTransform[0];
        rightMinLon=-180.0;
    }

    int left2 = right1,right2 = nXSize;

    auto pDriver=poDataset->GetDriver();

    //第一部分
    if(right1>0)
    {
        adfGeoTransform[0]=leftMinLon;
        auto poOutput1 = (GDALDataset *) GDALCreate(pDriver, outputFile1.c_str(), right1 - left1, nYSize,bandCount, dataType, nullptr);

        for (int i = 0; i < bandCount; i++)
        {
            GDALRasterBand *poBand = poDataset->GetRasterBand(i + 1);
            GDALRasterBand *poOutputBand = poOutput1->GetRasterBand(i + 1);

            // 分配内存来存储数据
            void *pData = malloc((right1 - left1) * nYSize * GDALGetDataTypeSize(dataType));
            if (pData == nullptr) {
                std::cerr << "Error: Failed to allocate memory for band data." << std::endl;
                GDALClose(poDataset);
                return false;
            }

            // 读取数据到内存
            poBand->RasterIO(GF_Read, left1, 0, right1 - left1, nYSize, pData, right1 - left1, nYSize,dataType, 0, 0);

            // 将数据写入到输出波段
            poOutputBand->RasterIO(GF_Write, 0, 0, right1 - left1, nYSize, pData, right1 - left1, nYSize,dataType, 0, 0);

            poOutputBand->SetNoDataValue(noDataValue);
            // 释放分配的内存
            free(pData);
        }

        poOutput1->SetGeoTransform(adfGeoTransform);
        GDALClose(poOutput1);
    }
    //第2部分
    if(right2>left2)
    {
        adfGeoTransform[0]=rightMinLon;
        auto poOutput2 = (GDALDataset *) GDALCreate(pDriver, outputFile2.c_str(), right2 - left2, nYSize, bandCount,dataType, nullptr);
        // 读取并保存第二个部分
        for (int i = 0; i <bandCount; i++)
        {
            GDALRasterBand *poBand = poDataset->GetRasterBand(i + 1);
            GDALRasterBand *poOutputBand = poOutput2->GetRasterBand(i + 1);

            // 分配内存来存储数据
            void *pData = malloc((right2 - left2) * nYSize * GDALGetDataTypeSize(dataType));
            if (pData == nullptr) {
                std::cerr << "Error: Failed to allocate memory for band data." << std::endl;
                GDALClose(poDataset);
                return false;
            }

            // 读取数据到内存
            poBand->RasterIO(GF_Read, left2, 0, right2 - left2, nYSize, pData, right2 - left2, nYSize,dataType, 0, 0);

            // 将数据写入到输出波段
            poOutputBand->RasterIO(GF_Write, 0, 0, right2 - left2, nYSize, pData, right2 - left2, nYSize,dataType, 0, 0);
            poOutputBand->SetNoDataValue(noDataValue);
            // 释放分配的内存
            free(pData);
        }
        poOutput2->SetGeoTransform(adfGeoTransform);
        GDALClose(poOutput2);
    }

    // 关闭输入数据集
    GDALClose(poDataset);
    // 删除输入文件
    if (std::remove(fileName) != 0)
    {
        std::cerr << "Error: Failed to delete input file " << fileName << std::endl;
        return false;
    }
    return true;
}

bool loadRaster(const char* fileName, double* &geoTransform, int &xSize, int &ySize,GDALDataType &dataType, void** &pDataList,double &noDataValue,int &bandCount)
{
    GDALAllRegister();

    //注册所有的驱动
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    GDALAllRegister();

    //GDAL的driver
    auto poSrcDS = (GDALDataset*)GDALOpen(fileName, GA_ReadOnly);
    if (poSrcDS == nullptr)
        return false;

    geoTransform= new double[6];
    poSrcDS->GetGeoTransform(geoTransform);

    xSize = poSrcDS->GetRasterXSize();
    ySize = poSrcDS->GetRasterYSize();

    dataType = poSrcDS->GetRasterBand(1)->GetRasterDataType();

    bandCount = poSrcDS->GetRasterCount();
    pDataList=new void*[bandCount];

    for (int bandId = 0; bandId < bandCount; bandId++)
    {
        // 分配数据空间
        size_t dataSize = GDALGetDataTypeSizeBytes(dataType) * xSize * ySize;
        pDataList[bandId]=CPLMalloc(dataSize);
        noDataValue=poSrcDS->GetRasterBand(bandId + 1)->GetNoDataValue();
        poSrcDS->GetRasterBand(bandId + 1)->RasterIO(GF_Read, 0, 0, xSize, ySize, pDataList[bandId],  xSize, ySize, dataType, 0, 0);

        if(!minMaxFlag)
        {
            double minMax[2];
            poSrcDS->GetRasterBand(bandId + 1)->ComputeRasterMinMax(0, minMax);  // 计算最小值和最大值
            minValue = minValue<minMax[0]?minValue:minMax[0];
            maxValue = maxValue>minMax[1]?maxValue:minMax[1];
        }
    }

    GDALClose(poSrcDS);

    return true;
}

#pragma endregion
/*********************************************tile相关操作***************************************************/
#pragma region tile与bound相关操作实现

#pragma region bound获取瓦片
__global__ void getGridTileCodeFromGridCode(const double* bound, double differ, int xSize, int ySize, int xStep, int yStep, int level, int originLevel, DGGSGridTileCode* tiffCodes, ICuDGGSKernel *pKernel)
{
    for ( size_t i = blockDim.x*blockIdx.x+threadIdx.x; i < xSize*ySize; i += blockDim.x*gridDim.x)
    {
        int x = int( i % xSize);
        int y = int( i / xSize);

        // 检查边界，防止越界访问
        if (x < xSize && y < ySize)
        {
            double lon = bound[2] + (x * xStep) * differ;
            double lat = bound[0] - (y * yStep) * differ;
            if (lon < -180) lon = 360.0 + lon;
            GeoCoord coord;
            coord.lon=lon;
            coord.lat=lat;
            DGGSElementCode tempCode;
            tempCode.level=level;
            tempCode.elementType=DGGSElementType::Cell;
            pKernel->encode(coord,tempCode);
            DGGSGridTileCode tileCode;
            gridCode2GridTiffCode(tempCode, tileCode, originLevel);
            tiffCodes[i]=tileCode;
        }
    }
}

bool getGridTileFromBound(const double *bound, int level, int originLevel, DGGSGridTileCode *&tileCodes, int &tiffCount) {
    //存储GridTIFF编码
    std::vector<DGGSGridTileCode> resultVector;

    double differ= calculateDiffer(level-1);

    int xSize=0;
    if(bound[3]-bound[2]>0)
        xSize=std::ceil((bound[3]-bound[2])/differ);
    else
        xSize=std::ceil((bound[2]-bound[3])/differ);
    int ySize=std::ceil((bound[0]-bound[1])/differ);


    int xStep=(int)(360.0/(pow(2,level-originLevel))/differ/10.0);
    int yStep=(int)(180.0/(pow(2,level-originLevel))/differ/10.0);


    double* bound_d= nullptr;
    cudaMalloc((void**)&bound_d,sizeof(double)*4);
    cudaMemcpy(bound_d,bound,sizeof(double)*4,cudaMemcpyHostToDevice);


    DGGSGridTileCode* gridTiffCodes_d= nullptr;


    // GPU 启动配置
    int testXSize=(int)std::ceil(xSize/xStep)+1;
    int testYSize=(int)std::ceil(ySize/yStep)+1;
    cudaMalloc((void**)&gridTiffCodes_d, sizeof(DGGSGridTileCode) * testXSize * testYSize);
    tileCodes=new DGGSGridTileCode[testXSize * testYSize];

    getGridTileCodeFromGridCode<<<pGridDimX, pBlockDimX>>>(bound_d, differ, testXSize, testYSize, xStep, yStep, level,
                                                           originLevel, gridTiffCodes_d, pCuKernel);
    cudaDeviceSynchronize();
    cudaMemcpy(tileCodes, gridTiffCodes_d, sizeof(DGGSGridTileCode) * testXSize * testYSize, cudaMemcpyDeviceToHost);

    cudaFree(bound_d);
    cudaFree(gridTiffCodes_d);

    for(int tileIndex=0;tileIndex<testXSize*testYSize;tileIndex++)
    {

        if(tileCodes[tileIndex].basePartition==0||tileCodes[tileIndex].basePartition==11)
            continue;
        int i=0;
        for(;i<resultVector.size();i++)
            if(tileCodes[tileIndex] == resultVector[i])
                break;
        if(i==resultVector.size())
            resultVector.push_back(tileCodes[tileIndex]);
    }
    delete[]tileCodes;


    //拷贝结果
    tiffCount = (int)resultVector.size();
    tileCodes = new DGGSGridTileCode[tiffCount];
    for (int i = 0; i < tiffCount; i++)
        tileCodes[i] = resultVector[i];

    return true;
}
#pragma endregion

#pragma region 瓦片获取bound

__global__ void code2Coord(DGGSElementCode * codes, int codeCount, GeoCoord* coords, ICuDGGSKernel*pKernel)
{
    for ( size_t i = blockDim.x*blockIdx.x+threadIdx.x; i < codeCount; i += blockDim.x*gridDim.x)
    {
        pKernel->decode(codes[i], coords[i]);
        if(coords[i].lon>=180)
            coords[i].lon-=360.0;
    }
}

DGGridDiamond getSpaceBoundFromGridTileCode(DGGSGridTileCode gridTiffCode, int originLevel)
{
    DGGridDiamond diamond;
    DGGSElementCode code;
    code.elementType=pElementType;

    gridTiffCode2GridCode(gridTiffCode, code, originLevel);
    diamond.basePartition=gridTiffCode.basePartition;

    auto rowMin = code.row;
    auto colMin = code.col;
    long long int rowMax;
    long long int colMax;
    if (code.level<originLevel)
    {
        rowMax = (long long)std::pow(2,code.level);
        colMax = (long long) std::pow(2,code.level);
    }
    else
    {
        rowMax = rowMin + (long long)std::pow(2, originLevel);
        colMax = colMin + (long long)std::pow(2, originLevel);
    }


    auto codes=new DGGSElementCode[4];
    for(int i=0;i<4;i++)
    {
        codes[i].basePartition=gridTiffCode.basePartition;
        codes[i].level=code.level;
        codes[i].elementType=DGGSElementType::Cell;
    }

    codes[0].row=rowMin;
    codes[0].col= colMax - 1;

    codes[1].row= rowMax - 1;
    codes[1].col= colMin;

    codes[2].row= rowMin;
    codes[2].col= colMin;

    codes[3].row= rowMax - 1;
    codes[3].col= colMax - 1;

    auto coords= new GeoCoord[4];

    DGGSElementCode*codes_d= nullptr;
    cudaMalloc((DGGSElementCode**)&codes_d, sizeof(DGGSElementCode) * 4);
    cudaMemcpy(codes_d,codes, sizeof(DGGSElementCode) * 4, cudaMemcpyHostToDevice);

    GeoCoord* coords_d= nullptr;
    cudaMalloc((GeoCoord**)&coords_d, sizeof(GeoCoord) * 4);

    int codeCount=4;
    code2Coord<<<1, 4>>>(codes_d, codeCount, coords_d, pCuKernel);
    cudaDeviceSynchronize();
    auto error=cudaMemcpy(coords,coords_d, sizeof(GeoCoord) * 4, cudaMemcpyDeviceToHost);
    if(error!=cudaSuccess)
        printf("cudaMemcpy error:%s\n",cudaGetErrorString(error));

    //上
    diamond.latMax = coords[0].lat;

    //下
    diamond.latMin = coords[1].lat;

    //左
    diamond.lonMin = coords[2].lon;

    //右
    diamond.lonMax = coords[3].lon;

    if(codes[0].basePartition<6)
        if(codes[0].row == 0)
        {
            diamond.lonMin= (codes[0].basePartition-1)*baseDiamondWidth+(-180.0)-36.0;
        }
    if(codes[0].basePartition>5)
        if (codes[1].col == 0)
        {
            diamond.lonMin = (codes[0].basePartition - 6) * baseDiamondWidth + (-144.0) - 36.0;
        }

    delete[]coords;
    delete[]codes;
    cudaFree(coords_d);
    cudaFree(codes_d);

    return diamond;
}

#pragma endregion

#pragma endregion
/*********************************************Raster2Grid***************************************************/
#pragma region Raster2Grid实现

DGGSDataType createMap(std::string* srcFileNames,int fileCount,int &level,int originLevel)
{
    raster2TileMap=new std::multimap<RasterIndex,DGGSGridTileCode>[1];

    auto bound=new double[4];
    DGGSDataType dataType;
    for(int i=0;i<fileCount;i++)
    {
        //读取raster
        auto geoTransform=new double[6];
        int xSize, ySize;
        GDALDataType gdalDataType;
        double noDataValue;
        int bandCount;
        void** pDataList= nullptr;
        loadRaster(srcFileNames[i].c_str(), geoTransform, xSize, ySize, gdalDataType, pDataList, noDataValue, bandCount);

        //获取Bound和层级
        getBoundFromGeoTrans(geoTransform, xSize, ySize,bound);
        //level=calculateLevel(geoTransform[1]);

        if(i==0)
            dataType=static_cast<DGGSDataType>(gdalDataType);

        //获取GridTIFF
        DGGSGridTileCode *gridTiffCodes= nullptr;
        int gridTiffCount=0;
        getGridTileFromBound(bound, level, originLevel, gridTiffCodes, gridTiffCount);

        for(int gridTiffIndex=0;gridTiffIndex<gridTiffCount;gridTiffIndex++)
            raster2TileMap->insert({srcFileNames[i], gridTiffCodes[gridTiffIndex]});

        delete[]gridTiffCodes;
        gridTiffCodes= nullptr;
    }
    delete[]bound;
    bound= nullptr;

    minMaxFlag= true;
    return dataType;
}


template<class T>
__device__ T getValueFromCoord(GeoCoord coord, const double* trans, int xSize, int ySize, T* srcValueList,double noValue)
{
    double lon=coord.lon;
    double lat=coord.lat;

    double temp = trans[1] * trans[5] - trans[2] * trans[4];
    double tempX = 0.0;
    double tempY = 0.0;
    tempX = (trans[5] * (lon - trans[0]) - trans[2] * (lat - trans[3])) / temp;
    tempY = (trans[1] * (lat - trans[3]) - trans[4] * (lon - trans[0])) / temp;

    //会把插值点的大致行列号拟合为所在格元左下角
    int x = static_cast<int>(tempX);
    int y = static_cast<int>(tempY);

    if (lon < trans[0] + trans[1] * xSize && x == xSize)
        x--;
    if (lat > trans[3] + trans[5] * ySize && y == ySize)
        y--;

    T pixValue;
    if (x >= xSize || y >= ySize || x < 0 || y < 0)
        pixValue = noValue;
    else
        pixValue = srcValueList[(int)(y * xSize + x)];
    return pixValue;
}


template<class T>
__global__ void setGridTIFFValuesFromRaster(T* srcValueList, double* trans, int xSize, int ySize, T* tagValueList, GridTransform tagTrans, int gridSize, ICuDGGSKernel* pKernel, double noValue)
{
    //格网编码转经纬度
    for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < gridSize * gridSize; i += blockDim.x * gridDim.x)
    {
        int row = (int)(i / gridSize);
        int col = (int)(i- row * gridSize);

        //获取属性值
        GeoCoord coord;

        DGGSElementCode code;
        code.elementType= DGGSElementType::Cell;
        code.level= tagTrans.level;
        code.row=row+tagTrans.rowOrigin;
        code.col=col+tagTrans.colOrigin;
        code.basePartition= tagTrans.basePartition;


        pKernel->decode(code,coord);

        T pixValue = getValueFromCoord(coord, trans, xSize, ySize, srcValueList,noValue);

        int x, y;
        code2XY(row, col, 0, 0, gridSize, x, y);

        if(tagValueList[(gridSize * 2 - 2 - y) * (gridSize * 2) + x] !=pixValue)
        {
            tagValueList[(gridSize * 2 - 2 - y) * (gridSize * 2) + x] = pixValue;
            tagValueList[(gridSize * 2 - 2 - y) * (gridSize * 2) + x + 1] = pixValue;
        }
    }
}


template<class T>
void singleRaster2Grid(const std::string& srcFileName, const char *exFormat,int originLevel,const std::string& tagFilePath)
{
    //根据raster2TileMap获取瓦片
    auto range = raster2TileMap->equal_range(srcFileName);
    std::vector<DGGSGridTileCode> gridTiffCodes;
    for (auto it = range.first; it != range.second; ++it)
        gridTiffCodes.push_back(it->second);

    //加载raster数据
    auto geoTransform=new double[6];
    int xSize, ySize;
    GDALDataType dataType;
    double noDataValue;
    int bandCount;
    void** pDataList= nullptr;
    loadRaster(srcFileName.c_str(), geoTransform, xSize, ySize, dataType, pDataList, noDataValue, bandCount);

    //修正nodata
    if(noDataValue>std::numeric_limits<T>::max()||noDataValue<std::numeric_limits<T>::min())
    {
        if(minValue>std::numeric_limits<T>::min())
            noDataValue=std::numeric_limits<T>::min();
        else if(maxValue<std::numeric_limits<T>::max())
            noDataValue=std::numeric_limits<T>::max();
    }//noDataValue=static_cast<T>(noDataValue);

    //GPU端raster数据
    T* srcValueList_d = nullptr;
    cudaMalloc((T**)&srcValueList_d, sizeof(T) * xSize *ySize);

    //坐标转换信息
    double* transform_d = nullptr;
    cudaMalloc((double**)&transform_d, sizeof(double) * 6);
    cudaMemcpy(transform_d, geoTransform, sizeof(double) * 6, cudaMemcpyHostToDevice);

    int gridSize=(int)pow(2,originLevel);
    auto dggsDataType=static_cast<DGGSDataType>(dataType);


    //遍历每个瓦片赋值（没有就创建加入）
    for(auto currentGridTiffCode:gridTiffCodes)
    {
        std::string gridTiffFileName=tagFilePath+"/"+getGridTiffString(currentGridTiffCode,originLevel)+exFormat;

        //设置格网转换参数
        DGGSElementCode code;
        gridTiffCode2GridCode(currentGridTiffCode,code,originLevel);
        GridTransform gridTransform{};
        gridTransform.level=code.level;
        gridTransform.rowOrigin=code.row;
        gridTransform.colOrigin=code.col;
        gridTransform.basePartition=code.basePartition;

        //设置GridTIFF
        GridTIFF gridTiff;

        //初始化GridTIFF
        if (!(std::filesystem::exists(gridTiffFileName)))
        {
            gridTiff.setSize(gridSize);
            gridTiff.setDataType(dggsDataType);
            gridTiff.setBandCount(bandCount);
            gridTiff.setTrans(gridTransform);
            gridTiff.setNoValue(noDataValue);
            gridTiff.setElementType(pElementType);
            gridTiff.setGridType(pGridType);
        }

        //获取GridTIFF数据
        else
        {
            gridTiff.loadFromFile(gridTiffFileName.c_str());
        }


        //GPU端GridTIFF数据
        T* tagValueLists_d = nullptr;
        cudaMalloc((T**)&tagValueLists_d, sizeof(T) * gridSize * gridSize*4);

        //逐个波段赋值
        for(int bandIndex=0;bandIndex<gridTiff.getBandCount();bandIndex++)
        {
            T *tagValueList=reinterpret_cast<T*>(gridTiff.getBandValue(bandIndex));
            T *srcValueList=reinterpret_cast<T*>(pDataList[bandIndex]);
            std::fill(tagValueList,tagValueList+gridSize * gridSize * 4,noDataValue);

            //CPU数据复制到GPU
            cudaMemcpy(tagValueLists_d, tagValueList, sizeof(T) * gridSize * gridSize * 4, cudaMemcpyHostToDevice);

            cudaMemcpy(srcValueList_d, srcValueList, sizeof(T) * xSize * ySize, cudaMemcpyHostToDevice);

            //GPU运算
            setGridTIFFValuesFromRaster<<<pGridDimX,pBlockDimX>>>(srcValueList_d, transform_d, xSize, ySize, tagValueLists_d, gridTransform, gridSize, pCuKernel, noDataValue);
            cudaDeviceSynchronize();

            //GPU数据复制到CPU
            cudaMemcpy(tagValueList, tagValueLists_d, sizeof(T) *  gridSize * gridSize * 4, cudaMemcpyDeviceToHost);
        }
        cudaFree(tagValueLists_d);

        gridTiff.saveToFile(gridTiffFileName.c_str());
    }

    cudaFree(srcValueList_d);
    cudaFree(transform_d);
    delete[]geoTransform;
    for(int i=0;i<bandCount;i++)
        CPLFree(pDataList[i]);
}


void raster2Grid(const char* srcFilePath,  const char* tagFilePath, const char*exFormat,int gridDimX, int blockDimX, DGGSGridType gridType, DGGSElementType elementType, int level, int originLevel)
{
    pGridDimX=gridDimX;
    pBlockDimX=blockDimX;

    //获取所有文件
    std::string * rasterFiles= nullptr;
    int rasterFileCount=0;
    if(fs::is_directory(srcFilePath))
        getAllFilesFromPath(srcFilePath,exFormat,rasterFiles,rasterFileCount);
    else if (fs::is_regular_file(srcFilePath)) {
        rasterFileCount=1;
        rasterFiles=new std::string[1];
        rasterFiles[0]=srcFilePath;
    }

    pGridType=gridType;
    pElementType=elementType;
    pCuKernel=CreateCuKernelObj(gridType);

    //创建索引表
    DGGSDataType dataType=createMap(rasterFiles,rasterFileCount,level,originLevel);


    #define SING_HANDLE_DATA_TYPE(T,D) \
        case D:                        \
        {                              \
            singleRaster2Grid<T>(rasterFiles[i],exFormat,originLevel,tagFilePath);   \
            break;                                  \
        }

    //singleRaster2Grid为瓦片赋值
    for(int i=0;i<rasterFileCount;i++)
    {
        switch (dataType) {
            SING_HANDLE_DATA_TYPE(unsigned char, DGGS_Byte)
            SING_HANDLE_DATA_TYPE(unsigned short, DGGS_UInt16)
            SING_HANDLE_DATA_TYPE(short, DGGS_Int16)
            SING_HANDLE_DATA_TYPE(unsigned int, DGGS_UInt32)
            SING_HANDLE_DATA_TYPE(int, DGGS_Int32)
            SING_HANDLE_DATA_TYPE(float, DGGS_Float32)
            SING_HANDLE_DATA_TYPE(double, DGGS_Float64)
            default:
                break;
        }
    }

    DestroyCuKernelObj(pCuKernel);
    cudaDeviceReset();
    raster2TileMap->clear();
    delete[]raster2TileMap;
    delete[]rasterFiles;
}


#pragma endregion
/***********************************************************************************************************/
#pragma region Grid2Raster实现

template<class T>
__global__ void setRasterValuesFromGrid(T* srcValueList, GridTransform srcTrans, int tileSize, T* tagValueList, double* trans, int xSize, int ySize, double noValue,double differTol,  ICuDGGSKernel* pKernel)
{
    for ( int i = blockDim.x*blockIdx.x+threadIdx.x; i < xSize*ySize; i += blockDim.x*gridDim.x)
    {
        int y=(int)(i/xSize);
        int x=(int)(i-y*xSize);

        if(x < 0 || x >= xSize || y < 0 || y > ySize - 1)
            continue;

        //geoTiff的经纬度坐标
        double lon = trans[0] + (x+0.5) * trans[1] + (y+0.5) * trans[2];
        double lat = trans[3] + (x+0.5) * trans[4] + (y+0.5) * trans[5];
        if(lon>180)
            lon=lon-differTol-360;
        if(lon<-180)
            lon=lon+differTol+360;

        DGGSElementCode code;
        code.elementType=DGGSElementType::Cell;
        code.level=srcTrans.level;
        GeoCoord coord;
        coord.lon=lon;
        coord.lat=lat;

        pKernel->encode(coord,code);
        //属性值
        if(srcTrans.basePartition!=code.basePartition||
           code.row<srcTrans.rowOrigin||code.row>=srcTrans.rowOrigin+tileSize||
           code.col<srcTrans.colOrigin||code.col>=srcTrans.colOrigin+tileSize
                )
        {
            tagValueList[y * xSize + x]=(T)noValue;
        }

        else
        {
            int gridX, gridY;
            code2XY(code.row, code.col, srcTrans.rowOrigin, srcTrans.colOrigin, tileSize, gridX, gridY);
            int width = tileSize * 2;
            int height = tileSize * 2;

            if (gridX >= 0 && gridX < width && gridY >= 0 && gridY < height - 1)
                tagValueList[y * xSize + x]=srcValueList[(height - 2 - gridY) * width + gridX];
            else
                tagValueList[y * xSize + x]=(T)noValue;
        }
    }
}

template<class T>
void singleGridTranslate(const char* srcFilePath,const char* gridTiffName,const char*exFormat, const char* tagFilePath)
{
    //raster文件
    std::string tagFileName=std::string(tagFilePath)+"/"+std::string(gridTiffName)+exFormat;

    //读取GridTIFF数据
    GridTIFF gridTiff;
    std::string gridTiffFileName=std::string(srcFilePath)+"/"+std::string(gridTiffName)+exFormat;
    if (!(std::filesystem::exists(gridTiffFileName)))
        return;
    gridTiff.loadFromFile(gridTiffFileName.c_str());

    int gridSize=gridTiff.getSize();

    //获取GridTIFF的转换参数
    GridTransform srcTrans=gridTiff.getTrans();
    T* srcValueList_d = nullptr;
    cudaMalloc((T**)&srcValueList_d, sizeof(T) * gridTiff.getSize() * gridTiff.getSize() * 4);

    int originLevel=(int)log2(gridTiff.getSize());

    //获取GridTIFFCode
    DGGSGridTileCode gridTiffCode=getGridTiffCode(gridTiffName);
    DGGridDiamond diamond=getSpaceBoundFromGridTileCode(gridTiffCode, originLevel);

    double differ= calculateDiffer(srcTrans.level-1);


    /*-180、180交界*/
    if( diamond.lonMin>diamond.lonMax&&abs(diamond.lonMin-diamond.lonMax)>72.0&&abs(diamond.lonMin-180)< differ)
        diamond.lonMin=diamond.lonMin-360.0;

    //可变范围
    double rangeDiffer=(diamond.lonMax-diamond.lonMin)/(gridSize*2-1);
    if(rangeDiffer<0.0)
        rangeDiffer=(diamond.lonMax-diamond.lonMin+360.0)/(gridSize*2-1)+2*differ;
    else
        rangeDiffer=(diamond.lonMax-diamond.lonMin)/(gridSize*2-1)+2*differ;
    diamond.lonMin=diamond.lonMin-rangeDiffer;
    diamond.lonMax=diamond.lonMax+rangeDiffer;
    diamond.latMax=diamond.latMax+rangeDiffer>90.0?90.0:diamond.latMax+rangeDiffer;
    diamond.latMin=diamond.latMin-rangeDiffer<-90.0?-90.0:diamond.latMin-rangeDiffer;

    //处理边界问题
    if(diamond.lonMin>diamond.lonMax)
    {
        diamond.lonMax=diamond.lonMax+360.0;
    }
    diamond.lonMin=std::ceil((diamond.lonMin-(-180.0))/differ)*differ+(-180.0);
    diamond.latMax=std::ceil((diamond.latMax-(90.0))/-differ)*(-differ)+(90.0);

    double geoTransform[6]={
            diamond.lonMin,differ,0,
            diamond.latMax,0,-differ};

    double bound[4]={diamond.latMax,diamond.latMin,diamond.lonMin,diamond.lonMax};

    if(bound[0]>mergeBound[0])
        mergeBound[0]=bound[0];
    if(bound[1]<mergeBound[1])
        mergeBound[1]=bound[1];
    if(bound[2]<mergeBound[2])
        mergeBound[2]=bound[2];
    if(bound[3]>mergeBound[3])
        mergeBound[3]=bound[3];
    mergeDiffer=differ;

    int tagYSize = std::round(((diamond.latMax - diamond.latMin) / differ));
    int tagXSize = std::round(((diamond.lonMax - diamond.lonMin) / differ));
    if(tagXSize<=0)
        tagXSize = std::round(((diamond.lonMax + 360.0 - diamond.lonMin) / differ));

    double differTol=0.0;
    if(diamond.lonMax>180.0)
    {
        differTol=-180.0+std::ceil(360.0/ differ)*differ-180.0;
    }

    int bandCount=gridTiff.getBandCount();
    double noValue=gridTiff.getNoValue();
    auto dataType=static_cast<GDALDataType>(gridTiff.getDataType());

    T **rasterValues = new T*[bandCount];

    //GPU端目标数据
    double* transform_d = nullptr;
    cudaMalloc((double**)&transform_d, sizeof(double) * 6);
    cudaMemcpy(transform_d, geoTransform, sizeof(double) * 6, cudaMemcpyHostToDevice);

    T* tagValueLists_d = nullptr;
    cudaMalloc((T**)&tagValueLists_d, sizeof(T) * tagYSize * tagXSize);

    for(int bandIndex=0;bandIndex<gridTiff.getBandCount();bandIndex++)
    {
        //CPU数据复制到GPU
        T *srcValueList=reinterpret_cast<T*>(gridTiff.getBandValue(bandIndex));

        rasterValues[bandIndex]=new T[tagXSize * tagYSize];
        std::fill(rasterValues[bandIndex], rasterValues[bandIndex] + tagXSize * tagYSize, noValue);
        cudaMemcpy(tagValueLists_d, rasterValues[bandIndex], tagXSize * tagYSize* sizeof(T), cudaMemcpyHostToDevice);

        cudaMemcpy(srcValueList_d, srcValueList, sizeof(T) * gridSize * gridSize * 4, cudaMemcpyHostToDevice);

        setRasterValuesFromGrid<<<pGridDimX,pBlockDimX>>>(srcValueList_d, srcTrans, gridTiff.getSize(), tagValueLists_d, transform_d, tagXSize, tagYSize, noValue,differTol,pCuKernel);
        auto error=cudaDeviceSynchronize();
        if (error != cudaSuccess)
        {
            std::cout << "setRasterValuesFromGrid error: " << error<<std::endl;
        }

        cudaMemcpy(rasterValues[bandIndex], tagValueLists_d, tagXSize * tagYSize* sizeof(T), cudaMemcpyDeviceToHost);
    }

    saveRaster(tagFileName.c_str(), geoTransform, tagXSize, tagYSize, dataType, rasterValues,noValue,bandCount);

    if(geoTransform[0]+tagXSize*geoTransform[1]>180.0||geoTransform[0]<-180.0)
    {
        clip(tagFileName.c_str());
    }

    cudaFree(srcValueList_d);
    cudaFree(transform_d);
    cudaFree(tagValueLists_d);
    for(int bandIndex=0;bandIndex<gridTiff.getBandCount();bandIndex++)
        delete []rasterValues[bandIndex];
    delete []rasterValues;
    cudaFree(0);
}

template<class T>
__global__ void mergeRasterValues(T* srcValueList, int srcXSize,int srcYSize,T*tagValueList,double noValue)
{
    for ( size_t i = blockDim.x*blockIdx.x+threadIdx.x; i < srcXSize*srcYSize; i += blockDim.x*gridDim.x)
    {
        if(srcValueList[i]!=noValue&&tagValueList[i]==noValue)
        {
            tagValueList[i]=srcValueList[i];
        }
    }
}

template<class T>
void mergeRaster(const char* srcFilePath,const char* tagFileName,const char*exFormat)
{
    std::string *rasterFileNames= nullptr;
    int rasterTiffCount=0;
    getAllFilesFromPath(srcFilePath,exFormat, rasterFileNames, rasterTiffCount);

    //合并图层
    if(mergeBound[3]>180.0)mergeBound[3]=180.0;
    if(mergeBound[2]<-180.0)mergeBound[2]=-180.0;
    int mergeXSize = std::ceil((mergeBound[3] - mergeBound[2]) / mergeDiffer);
    int mergeYSize = std::ceil((mergeBound[0] - mergeBound[1]) / mergeDiffer);
    double mergeTrans[6]={
            mergeBound[2],mergeDiffer,0,
            mergeBound[0],0,-mergeDiffer};

    //合并数据
    GDALDataset* mergeSrcDS = nullptr;
    //获取值
    for(int i=0; i < rasterTiffCount; i++)
    {
        //读取当前Raster数据
        auto currentTrans=new double[6];
        int currentXSize=0;
        int currentYSize=0;
        void** currentValues= nullptr;
        double noDataValue=0;

        std::cout<<i<<"\t"<<rasterFileNames[i]<<std::endl;
        loadRaster(rasterFileNames[i].c_str(),currentTrans, currentXSize,currentYSize,mergeDataType, currentValues,noDataValue,mergeBandCount);

        //初始化合并图层
        if(i==0)
        {
            T** mergeValues=new T*[mergeBandCount];
            for(int mergeBandId=0;mergeBandId<mergeBandCount;mergeBandId++)
            {
                mergeValues[mergeBandId]= new T[mergeXSize * mergeYSize];
                std::fill(mergeValues[mergeBandId], mergeValues[mergeBandId] + mergeXSize * mergeYSize, (T)noDataValue);
            }

            saveRaster(tagFileName, mergeTrans, mergeXSize, mergeYSize, mergeDataType, mergeValues,noDataValue,mergeBandCount);
            for (int mergeBandId = 0; mergeBandId < mergeBandCount; mergeBandId++)
                delete[]mergeValues[mergeBandId];
            delete[]mergeValues;
            mergeSrcDS= (GDALDataset*)GDALOpen(tagFileName, GA_Update);
        }

        //计算偏移
        int xOffset = (int)round((currentTrans[0] - mergeBound[2]) / mergeDiffer);
        int yOffset = (int)round((mergeBound[0] - currentTrans[3]) / mergeDiffer);

        //合并图层数据
        T*mergeValue_d= nullptr;
        cudaMalloc((T**)&mergeValue_d, sizeof(T) * currentXSize * currentYSize);

        T*currentValue_d= nullptr;
        cudaMalloc((T**)&currentValue_d, sizeof(T) * currentXSize * currentYSize);
        GDALDataType gdalDataType=mergeSrcDS->GetRasterBand(1)->GetRasterDataType();

        T* mergeValue=new T[currentXSize * currentYSize];

        //逐个波段赋值
        for(int bandId=0;bandId<mergeBandCount;bandId++)
        {
            mergeSrcDS->GetRasterBand(bandId + 1)->RasterIO(GF_Read, xOffset, yOffset,currentXSize,currentYSize, mergeValue,currentXSize,currentYSize, gdalDataType, 0, 0);

            //数据复制到GPU
            cudaMemcpy(currentValue_d,currentValues[bandId],sizeof(T) * currentXSize * currentYSize,cudaMemcpyHostToDevice);
            cudaMemcpy(mergeValue_d,mergeValue,sizeof(T) * currentXSize * currentYSize,cudaMemcpyHostToDevice);

            mergeRasterValues<<<pGridDimX,pBlockDimX>>>(currentValue_d,currentXSize,currentYSize,mergeValue_d,(T)noDataValue);
            cudaDeviceSynchronize();

            //数据复制到CPU
            cudaMemcpy(mergeValue, mergeValue_d, sizeof(T) * currentXSize * currentYSize, cudaMemcpyDeviceToHost);
            mergeSrcDS->GetRasterBand(bandId + 1)->RasterIO(GF_Write, xOffset, yOffset,currentXSize,currentYSize, mergeValue,currentXSize,currentYSize, gdalDataType, 0, 0);
        }

        cudaFree(currentValue_d);
        cudaFree(mergeValue_d);
        delete[]currentTrans;
        for(int bandIndex=0;bandIndex<mergeBandCount;bandIndex++)
            CPLFree(currentValues[bandIndex]);
        delete[]currentValues;
        delete[]mergeValue;
    }

    GDALClose((GDALDatasetH)mergeSrcDS);

    //删除源raster文件
    delete[]rasterFileNames;
}

void singleGrid2Raster(const char* srcFilePath, const char* tagFilePath,const char*exFormat)
{
    DGGSDataType dataType=DGGSDataType::DGGS_Unknown;
    int level=0;
    minMaxFlag=true;

    std::string *gridTiffNames= nullptr;
    int gridTiffCount=0;
    getAllFilesFromPath(srcFilePath,exFormat,gridTiffNames,gridTiffCount);
    if(gridTiffCount==0)
    {
        std::cerr<<"no GridTIFF exists"<<std::endl;
        return;
    }
    else
    {
        //获取GridTIFF参数
        GridTIFF gridTiff;
        gridTiff.loadFromFile(gridTiffNames[0].c_str());
        dataType=gridTiff.getDataType();
        level=gridTiff.getTrans().level;
        mergeDiffer= calculateDiffer(level-1);
        pGridType=gridTiff.getGridType();
        pElementType=gridTiff.getElementType();
        pCuKernel= CreateCuKernelObj(pGridType);
    }

    //获取文件名
    for(int i=0;i<gridTiffCount;i++)
    {
        size_t left = gridTiffNames[i].find_last_of("/\\");
        size_t right = gridTiffNames[i].find_last_of('.');
        gridTiffNames[i] = (right != std::string::npos) ? gridTiffNames[i].substr(left+1,right-left-1) : gridTiffNames[i];
    }

    std::string mergeFileName=tagFilePath;
    mergeFileName=mergeFileName+exFormat;

    //转换
    for(int i=0;i<gridTiffCount;i++)
    {

        std::cout<<i<<"\t"<<gridTiffNames[i]<<std::endl;


#define TRANS_HANDLE_DATA_TYPE(T, D) \
            case D:                        \
            {                          \
                gridTiffNames[i]=fs::path(gridTiffNames[i]).stem().string();        \
                singleGridTranslate<T>(srcFilePath, gridTiffNames[i].c_str(),exFormat, tagFilePath);   \
                break;                                  \
            }
        switch (dataType) {
            TRANS_HANDLE_DATA_TYPE(unsigned char, DGGS_Byte)
            TRANS_HANDLE_DATA_TYPE(unsigned short, DGGS_UInt16)
            TRANS_HANDLE_DATA_TYPE(short, DGGS_Int16)
            TRANS_HANDLE_DATA_TYPE(unsigned int, DGGS_UInt32)
            TRANS_HANDLE_DATA_TYPE(int, DGGS_Int32)
            TRANS_HANDLE_DATA_TYPE(float, DGGS_Float32)
            TRANS_HANDLE_DATA_TYPE(double, DGGS_Float64)
            default:
                break;
        }
    }

    //超过13级不进行合并
     if(level>MaxMergeLevel)
    {
        delete[]gridTiffNames;
        return;
    }

    //合并
#define MERGE_HANDLE_DATA_TYPE(T, D) \
        case D:                        \
        {                          \
            mergeRaster<T>(tagFilePath,mergeFileName.c_str(),exFormat); \
            break;                                  \
        }
    switch (dataType)
    {
        MERGE_HANDLE_DATA_TYPE(unsigned char, DGGS_Byte)
        MERGE_HANDLE_DATA_TYPE(unsigned short, DGGS_UInt16)
        MERGE_HANDLE_DATA_TYPE(short, DGGS_Int16)
        MERGE_HANDLE_DATA_TYPE(unsigned int, DGGS_UInt32)
        MERGE_HANDLE_DATA_TYPE(int, DGGS_Int32)
        MERGE_HANDLE_DATA_TYPE(float, DGGS_Float32)
        MERGE_HANDLE_DATA_TYPE(double, DGGS_Float64)
        default:
            break;
    }
    delete[]gridTiffNames;

    if (fs::exists(tagFilePath) && fs::is_directory(tagFilePath))
                fs::remove_all(tagFilePath);

    //还原变量
    delete raster2TileMap;
    mergeBound[0]=-90;
    mergeBound[1]=90;
    mergeBound[2]=180;
    mergeBound[3]=-180;
    mergeBandCount=0;
    mergeDiffer=0.0;
    mergeDataType=GDT_Unknown;
}

void multGrid2Raster(const char* srcFilePath, const char* tagFilePath,const char*exFormat)
{
    std::string* fileDirectories= nullptr;
    int fileDirectoriesCount=0;

    getAllDirectoriesFromPath(srcFilePath,fileDirectories,fileDirectoriesCount);

    for(int i=0;i<fileDirectoriesCount;i++)
    {
        std::string currentSrcDir=fileDirectories[i];
        fs::path p(currentSrcDir);
        std::string currentTagDir=std::string(tagFilePath)+"/"+p.filename().string();
        if (!fs::exists(currentTagDir))   // 检查目录是否存在
            fs::create_directories(currentTagDir);

        singleGrid2Raster(currentSrcDir.c_str(), currentTagDir.c_str(),exFormat);
    }
    delete[]fileDirectories;
}

void grid2Raster(const char* srcFilePath, const char* tagFilePath,const char*exFormat,int gridDimX,int blockDimX,int flag)
{
    pGridDimX=gridDimX;
    pBlockDimX=blockDimX;

    fs::create_directory(tagFilePath);

    if(1==flag)
        singleGrid2Raster(srcFilePath, tagFilePath,exFormat);
    else
        multGrid2Raster(srcFilePath,tagFilePath,exFormat);
    DestroyCuKernelObj(pCuKernel);
    cudaDeviceReset();
}

#pragma endregion